function [hessian,scoreQML,scoreSMM,dQ,dQMLpart,dSMMpart,errorMes] = getScoreHessian(paramsValues,setupFilter,outFilter,stepOption)

errorMes = 0;
MPon     = 0;
% Computing dPyy and dInov by numerical differentiation
numParams      = length(setupFilter.selectParams);
[samples,dimy] = size(setupFilter.data);
AA             = zeros(numParams,numParams);
hessian        = zeros(numParams,numParams);
scoreQML       = zeros(samples,numParams);
scoreSMM       = zeros(samples,numParams);
dQ             = zeros(numParams,1);
dQMLpart       = zeros(numParams,1);
dSMMpart       = zeros(numParams,1);
dPyy           = zeros(dimy,dimy,samples,numParams);
dInov          = zeros(dimy,samples,numParams);
epsValue       = setupFilter.epsValue;
if isfield(outFilter,'momCondition') == 1
    gradMom = zeros(length(outFilter.momCondition),numParams);
end

%epsValueDefault = epsValue;
if MPon == 1
    parfor i=1:numParams
        %epsValue = epsValueDefault*abs(paramsValues(i));
        [score_i,dPyy_i,dInov_i,gradMom_i] = getRowOfScore_mp(paramsValues,setupFilter,outFilter,epsValue,stepOption,i);
        dPyy(:,:,:,i)  = dPyy_i;
        dInov(:,:,i)   = dInov_i;
        gradMom(:,i)   = gradMom_i;
    end

    % Calculating the Hessian matrix
    parfor t=setupFilter.numInitial+1:samples
    %for t=setupFilter.numInitial+1:samples
       AA = getHessianPerPeriond_mp(setupFilter,outFilter,numParams,dPyy,dInov,t);
       hessian = hessian + AA;
    end
else
    for i=1:numParams
        i
        % Positive step
        paramsValues_epsp = paramsValues;
        paramsValues_epsp(i) = paramsValues_epsp(i) + epsValue;
        [Q_epsp,errorMes_epsp,model_epsp,~,outFilter_epsp] = objectFuncQML(paramsValues_epsp,setupFilter);
        if errorMes_epsp == 1
            disp(['No positive gradient step at position ',num2str(i)])
            errorMes = 1;
            return
        end
        
        % Negative step
        if stepOption == 2
            paramsValues_epsm = paramsValues;
            paramsValues_epsm(i) = paramsValues_epsm(i) - epsValue;
            [Q_epsm,errorMes_epsm,~,~,outFilter_epsm] = objectFuncQML(paramsValues_epsm,setupFilter);
            if errorMes_epsm == 1
                disp(['No negative gradient step at position ',num2str(i)])
                errorMes = 1;
                return
            end
        end
        
        if stepOption == 1
            scoreQML(setupFilter.numInitial+1:samples,i) = (outFilter_epsp.LogL(setupFilter.numInitial+1:end) - outFilter.LogL(setupFilter.numInitial+1:end))/epsValue;
            for t=1:samples
                if setupFilter.inclSurveys == 1
                    numObs = sum(~isnan(setupFilter.data(t,:)));
                elseif setupFilter.inclSurveys == 0
                    %numObs = sum(~isnan(setupFilter.data(t,1:end-length(model_epsp.matConShortRate))));
                    numObs  = sum(~isnan(outFilter_epsp.inov(:,t)));
                end
                dPyy(1:numObs,1:numObs,t,i) = (outFilter_epsp.Pyy(1:numObs,1:numObs,t) - outFilter.Pyy(1:numObs,1:numObs,t))/epsValue;
                dInov(1:numObs,t,i)         = (outFilter_epsp.inov(1:numObs,t)         - outFilter.inov(1:numObs,t))/epsValue;
            end
            if isfield(outFilter,'momCondition') == 1
                % The score of the moment conditions
                gradMom(:,i) = (outFilter_epsp.momCondition-outFilter.momCondition)/epsValue;
            end
            dQ(i,1)       = (Q_epsp-outFilter.Q)/epsValue;
            dQMLpart(i,1) = (outFilter_epsp.QMLpart-outFilter.QMLpart)/epsValue;
            dSMMpart(i,1) = (outFilter_epsp.SMMpart-outFilter.SMMpart)/epsValue;
        elseif stepOption == 2
            scoreQML(setupFilter.numInitial+1:samples,i) = (outFilter_epsp.LogL(setupFilter.numInitial+1:end) - outFilter_epsm.LogL(setupFilter.numInitial+1:end))/(2*epsValue);
            for t=1:samples
                if setupFilter.inclSurveys == 1
                    numObs = sum(~isnan(setupFilter.data(t,:)));
                elseif setupFilter.inclSurveys == 0
                    %numObs = sum(~isnan(setupFilter.data(t,1:end-length(model_epsp.matConShortRate))));
                    numObs  = sum(~isnan(outFilter_epsp.inov(:,t)));
                end
                dPyy(1:numObs,1:numObs,t,i) = (outFilter_epsp.Pyy(1:numObs,1:numObs,t) - outFilter_epsm.Pyy(1:numObs,1:numObs,t))/(2*epsValue);
                dInov(1:numObs,t,i)         = (outFilter_epsp.inov(1:numObs,t)         - outFilter_epsm.inov(1:numObs,t))/(2*epsValue);
            end
            if isfield(outFilter,'momCondition') == 1
                % The score of the moment conditions
                gradMom(:,i) = (outFilter_epsp.momCondition-outFilter_epsm.momCondition)/(2*epsValue);
            end
            dQ(i,1)       = (Q_epsp-Q_epsm)/(2*epsValue);
            dQMLpart(i,1) = (outFilter_epsp.QMLpart-outFilter_epsm.QMLpart)/(2*epsValue);
            dSMMpart(i,1) = (outFilter_epsp.SMMpart-outFilter_epsm.SMMpart)/(2*epsValue);
        end
   
    end
    % Calculating the Hessian matrix
    for t=setupFilter.numInitial+1:samples
        if setupFilter.inclSurveys == 1
            numObs = sum(~isnan(setupFilter.data(t,:)));
        elseif setupFilter.inclSurveys == 0
            %numObs = sum(~isnan(setupFilter.data(t,1:end-length(model_epsp.matConShortRate))));
            numObs  = sum(~isnan(outFilter_epsp.inov(:,t)));
        end
        for i=1:numParams
            for j=1:numParams
                tmp = outFilter.invPyy(1:numObs,1:numObs,t)*dPyy(1:numObs,1:numObs,t,i)*outFilter.invPyy(1:numObs,1:numObs,t)*dPyy(1:numObs,1:numObs,t,j);
                AA(i:i,j:j) = 0.5*trace(tmp)+dInov(1:numObs,t,i)'*outFilter.invPyy(1:numObs,1:numObs,t)*dInov(1:numObs,t,j);
            end
        end
        hessian = hessian + AA;
    end
end

if isfield(outFilter,'momCondition') == 1
    % Adding the moment conditions to the score 
    if setupFilter.SMMwithCSreg == 1 || setupFilter.SMMwithCSreg == 2
        scoreSMM = zeros(samples,numParams);
        for t=1:samples-setupFilter.CSregress_m
            scoreSMM(setupFilter.CSregress_m+t,:) = 2*setupFilter.lambdaSMM*(outFilter.modelMoments-setupFilter.dataMoments_mat(t,:))*setupFilter.W*gradMom;
        end
    else
        scoreSMM = zeros(samples,numParams);
        for t=1:samples
            scoreSMM(t,:) = 2*setupFilter.lambdaSMM*(outFilter.modelMoments-setupFilter.dataMoments_mat(t,:))*setupFilter.W*gradMom;
        end
        % Thus, scoreSMM/T = 2*setupFilter.lambdaSMM*outFilter.momCondition*setupFilter.W*gradMom        
    end
    % Adding the moment conditions to the Hessian
    hessian = hessian/(samples-setupFilter.numInitial) - 2*setupFilter.lambdaSMM*gradMom'*setupFilter.W*gradMom;
else
    scoreSMM = scoreQML*NaN;
    hessian  = hessian/(samples-setupFilter.numInitial);
end
if errorMes == 1
    scoreQML = scoreQML*NaN;
    scoreSMM = scoreSMM*NaN;
    hessian  = hessian*NaN;
end